
import sys
import indexing.featureindex


def printSnpEntries(details, sample):
  for base in details.cts.keys():
    if base != details.genomic and base != 'N':
      printSnpEntry(details, sample, base)
  
def printSnpEntry(details, sample, base):
  # Sample
  cols = [sample]
  # Unique Best?
  cols.append('Y')
  # Chromosome
  cols.append(details.chr.replace("chr", ""))
  # Position
  cols.append(details.spos)
  # Original Base
  cols.append(details.genomic)
  # Called Base
  cols.append(base)
  # Calling Reads
  cols.append(details.cts[base])

  # Parsing info string, of the form (1F23+1F30+1R20+1R21+1R25)
  readPositions = details.info[base][1:-1].split("+")
  startPositions = {}
  forwardCts = reverseCts = 0
  leftExtent = rightExtent = 0
  firstPosCts = lastThreePosCts = 0
  for pos in readPositions:
    brk = max(pos.find("R"), pos.find("F"))
    cts = int(pos[:brk])
    seqPos = int(pos[brk+1:])
    if pos[brk] == "F":
      startPositions[seqPos] = True
      forwardCts += cts
      leftExtent = max(leftExtent, seqPos)
      rightExtent = max(rightExtent, 31-seqPos)
    else:
      startPositions[31-seqPos] = True
      reverseCts += cts
      leftExtent = max(leftExtent, 31-seqPos)
      rightExtent = max(rightExtent, seqPos)
    if seqPos == 0: firstPosCts += cts
    if seqPos >= 29: lastThreePosCts += cts

  # Forward Calling Reads
  cols.append(forwardCts)
  # Reverse Calling Reads
  cols.append(reverseCts)
  # Calling Reads in Both Directions?
  if forwardCts>0 and reverseCts>0:
    cols.append('Y')
  else:
    cols.append('N')
  
  # Coverage
  cols.append(details.reads)
  # Fraction of Coverage Calling SNP
  cols.append("%.6f" % (1.0*details.cts[base]/details.reads))
  # Calling Reads Passing Filter
  cols.append(details.cts[base])
  # Fraction of Calling Reads Passing Filter
  cols.append("1.000000")
  # Coverage Passing Filter
  cols.append(details.reads)
  # Fraction of Coverage Passing Filter
  cols.append("1.000000")
  # Calling Read Start Positions
  cols.append(len(startPositions))
  # Extent of Left Calling Coverage
  cols.append(leftExtent)
  # Extent of Right Calling Coverage
  cols.append(rightExtent)
  # Extent of Calling Coverage
  cols.append(leftExtent + rightExtent + 1)
  # Fraction of Calling Reads with SNP in First Position
  cols.append("%.6f" % (1.0*firstPosCts/details.cts[base]))
  # Fraction of Calling Reads with SNP in Last Three Positions
  cols.append("%.6f" % (1.0*lastThreePosCts/details.cts[base]))

  print "\t".join([str(x) for x in cols])

def printUsage(args, usage, minArgs=2):
  if len(args) < minArgs or args[1].startswith("-h") or args[1].startswith("--h"):
    print usage
    sys.exit(1)

usageString = """
Usage: python %s COVERAGE_SUMMARY_FILE [CHR[:MIN-MAX]]
"""

if __name__=="__main__":
    printUsage(sys.argv, usageString % sys.argv[0], 2)
    # index = indexing.featureindex.FeatureIndex(sys.argv[1])
    coverageFile = sys.argv[1]
    (qchr, qmin, qmax) = (None, None, None)
    if len(sys.argv)>2:
      query = sys.argv[2]
      if query.find(":")>=0:
        qchr = query.split(":")[0]
        (qmin, qmax) = tuple(query.split(":")[1].split("-"))
        qmin = long(qmin)
        qmax = long(qmax)
      else:
        qchr = query
    query = (qchr, qmin, qmax)
    sample = coverageFile.split("/")[-1].split(".")[0]
    print >> sys.stderr, "# Coverage file =", coverageFile
    print >> sys.stderr, "# Sample =", sample
    print >> sys.stderr, "# Query =", query
    for feature in indexing.featureindex.getFeatures(coverageFile, query):
        details = feature.getDetails()
        if not details.snp: continue
        printSnpEntries(details, sample)

